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Abstract 

A class of high-order lowpass filters, the discrete singular convolution (DSC) 
filters, is utilized to facilitate the Fourier pseudospectral method for the so- 
lution of hyperbolic conservation law systems. The DSC filters are imple- 
mented directly in the Fourier domain (i.e., windowed Fourier pseudospec- 
tral method), while a physical domain algorithm is also given to enable the 
treatment of some special boundary conditions. By adjusting the effective 
wavenumber region of the DSC filter, the Gibbs oscillations can be removed 
effectively while the high resolution feature of the spectral method can be 
retained. The utility and effectiveness of the present approach is validated by 
extensive numerical experiments. 
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singular convolution filters; Gibbs oscillations 
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I. INTRODUCTION 



It is well known that due to their accuracy and efficiency, spectral methods have great 
advantages over local methods in solving applicable scientific and engineering problems 
[4,7,10,12,39]. Given a hyperbolic system of nonlinear conservation laws 



its solution may not exist in the classical sense because of possible discontinuities. The 
direct use of spectral methods to this problem will encounter the Gibbs oscillations [16], 
which lead to blow-up in the time integration. Therefore, it has been of great interest in 
making spectral methods applicable to the hyperbolic conservation law systems in the past 
two decades. The expectation, as point out by Gottlieb [15], is that the use of spectral 
methods will enable us not merely to capture the shock, but also to capture the delicate 
features and structures of the fiow. Obviously, spectral approaches will be extremely valuable 
to a class of aerodynamical problems that involve the interaction of both turbulence and 
shock [31]. Such an interaction is some of the most challenging problems in computation 
fiuid dynamics. 

The main objective of previous studies is to recover smooth solutions from those that are 
contaminated by Gibbs oscillations and, meanwhile, to improve the rate of convergence. For 
shock-capturing, many up-to-date local methods, such as weighted essentially non-oscillatory 
(WENO) scheme [29,35,36], central schemes [3,25,30,34], arbitrary-order non-oscillatory ad- 
vection scheme [38], Gas kinetic [47], and image processing based schemes [18,44], perform 
well. The success of these local shock- capturing schemes lies in their appropriate amount 
of intrinsic numerical dissipation, which is introduced either by explicit artificial viscosity 
terms, or by upwinding, or by appropriate local average (non-oscillatory central schemes). 




(1) 



with initial condition 



u{x, 0) = Uo{x), 



(2) 
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or by relaxation [24]. The characteristic decomposition based on Roe's mean matrix can 
also be considered as a local averaging on the Jacobian matrix A(u) = di{u)/du. In his 
recent work, LeVeque [28] outlined the relation between some approximate Riemann solvers 
and relaxation schemes. The relative amount of numerical dissipations may explain why 
the local characteristic decomposition is not necessary in low-order methods while it seems 
to be indispensable in high-order methods [35]. However, when these local methods are 
used in the cases where flow structures with flne details are needed to be resolved together 
with shocks, numerical dissipations are usually found to be so large that they smear the flne 
details [31]. The spectral methods, on the contrary, contribute very little numerical dissipa- 
tion and dispersion in principle when used to approximate spatial derivatives. Nevertheless, 
when they are applied to the approximation of spatial derivatives on a domain containing 
discontinuities, the Gibbs oscillations again have to be suppressed by appropriate means, for 
which there are two general routines: (1) explicit artiflcial viscosity, e.g. spectral viscosity 
(SV) method proposed by Tadmor [37] and (2) flltering, which is the central issue of the 
present study. 

At present, most fllters are constructed in the spectral domain and they can be called 
spectral fllters. Hussaini et al. [22] suggested a list of typical spectral fliers, i.e. Lanzos 
filter, raised cosine filter, sharpened raised cosine filer and exponential cutoff filter. More 
sophisticate and effective filters in spectral domain are Vandeven's p-th order filter [40] , Cai 
et aVs sawtooth function [5] and Gottlieb and Tadmor 's regularized Dirichlet function [13]. 

Apart from methods that are implemented in the spectral domain, filters in the physical 
domain are also developed. To design an appropriate filter in the physical domain is gener- 
ally more difficult than in the Fourier domain. A straightforward procedure is to make use 
of numerical dissipation partially contained in /ii^f/i-order shock-capturing schemes [6], such 
as the ENO scheme, where, actually, numerical dissipations were introduced both in the 
Fourier domain (via an exponential filter) and in the physical domain (via ENO polynomial 
interpolation). Such a strategy was generalized by Yee et al. [49] to the so-called charac- 
teristic filter method (with finite difference or finite volume) which has been successfully 
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applied to shock-capturing and turbulence simulations. Another procedure was initiated by 
Gottlieb [14] by using the Gegenbauer polynomial to resolve the oscillatory partial Fourier 
summation. Some promising numerical results from the filter approach can be found in Refs. 
[8,9,17,21]. 

Nevertheless, some prominent issues may hinder the success of filtered spectral methods 
for practical numerical computations involving shocks. In our view, there are two very 
important issues that are vita to the success of a filtered spectral scheme. The first issue is 
how to select optimal filters for shock-capturing. The issue is very complex and cannot have 
a unique answer at this point because there are so many properties to mind, such as flatness, 
ripple, filter length, effective frequency range and length of transition band, to name only 
a few. Loosely speaking, it is desirable to have filters that are free of dispersion errors, 
flat while having very small transition band, short in length while having high resolution. 
Moreover, what adds to the complexity is that solutions to different hyperbolic conservation 
law systems may have different Fourier spectral distributions. Consequently, one faces the 
difficulty that an optimal filter for a problem may not even work for another one. Therefore, 
it is desirable to have a filter which accounts for this change by an adjustable parameter. 

The second issue also concerns the change of solution, but is due to the time evolution 
for a given hyperbolic conservation law system. As the Fourier spectral distribution of a 
given problem changes with time, the question is how to implement the filter? It will be 
too dissipative for many systems if a filter is applied at each time step of the integration. 
Therefore, an adaptive implementation, which is controlled by a sensor during the time 
integration, is appropriate. Given the complexity, it is unlikely that there will be an ideal 
solution to these issues in the immediate future. Obviously, these challenging and interesting 
problems call for the further study of spectral filter approaches. 

The objective of this work is to construct, analyze and implement a class of high-order 
lowpass filters which can be used either in the Fourier domain or in the physical domain. 
These filters arc constructed within the framework of the discrete singular convolution (DSC) 
algorithm [41-43] and they are used to suppress Gibbs oscillations arising in the Fourier 
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pseudospectral method (FPM). Such lowpass filters have a wide effective wavenumber range 
that makes it possible to capture the fine flow structures, which is a desirable objective of 
spectral methods for hyperbolic conservation law systems. It is this effective wavenumber 
range that controls the resolution of the overall method. In our design, this effective range 
can be varied according to the resolution requirement by a parameter. Different DSC kernels 
have different magnitude responses and adjustability in the Fourier domain, which in turn 
influence the accuracy and resolution. The DSC filters used in this work were originally 
designed in a conjugate filter oscillation reduction (CFOR) [19,45,46] scheme. By 'conjugate 
filters' it means that the effective wavenumber range of the lowpass filter is largely overlapped 
with that of the highpass filter which is used for the approximation of spatial derivatives. 
The resulting Fourier domain algorithm is essentially a windowed Fourier pseudospectral 
method for shock-capturing. 

The rest of this paper is organized as follows. In Section II, a brief introduction will be 
devoted to the construction of the DSC lowpass filters followed by a discussion of two types 
of implementations of these filters to spectral methods. Extensive numerical experiments 
are carried out in Section III, including Burgers' equation, shock-tube problems, shock- 
entropy wave interaction, shock-vortex interaction. The shock-capturing ability and the 
high resolution feature of the present method are addressed. Concluding remarks end this 
paper. 



where T{x) is a singular kernel and rj^x) is an element of the space of test functions. Inter- 
esting singular kernels include that of Hilbert type, Abel type and delta type. The former 

5 



II. THEORY AND ALGORITHM 



A. DSC lowpass filter 



In the context of distribution theory, a singular convolution is defined by 




(3) 



two play important roles in the theory of analytical functions, processing of analytical sig- 
nals, theory of linear responses and Radon transform. Since delta type kernels are the key 
element in the theory of approximation and the numerical solution of differential equations, 
we focus on singular kernels of the delta type 

T{x) ^ 5^'^\x) , g = 0,1,2,- •• , (4) 

where superscript {q) denotes the gth-order 'derivative' of the delta distribution d{x), with 
respect to x, which should be understood as generalized derivatives of distributions. When 
q — 0, the kernel T{x) — 5{x) is important for the interpolation of functions. In this work, 
only the case of g = will be involved. One has to find appropriate approximations to 
the above singular kernel, which can not be directly realized in computers. To this end, we 
consider a sequence of approximations 

lim 5^f{x) = , g = 0, 1, 2, • • • (5) 

where a is a parameter which characterizes the approximation with the ckq being a gener- 
alized limit. Among various candidates of approximation kernels [42], regularized Shannon 
kernel (RSK) 

X ( X sinf(x-a:fc) / {x - Xkf\ 
'^-^^^ - "'^^ = l{x-x,) [ 2^ J 

is frequently used due to its close link to the Shannon sampling theorem. In this formula, A 

is the grid spacing and a determines the width of the Gaussian envelop. For a given cr 7^ 0, 

the limit of A — > reproduces the delta kernel (distribution). With the RSK, a function u 

can be approximated by a discrete convolution 

w 

u{x)^ J2 (^a,A(^-^k)u{xk) , q^0,l,2,--- , (7) 

k=-W 

where {xk}kL-w ^ of discrete grid points which are centered at x, and 21^ -|- 1 is 
the computational bandwidth, or effective kernel support, which is usually smaller than the 
computational bandwidth of the spectral method - the entire domain span. Equation (7) is 
often referred as a DSC and the RSK is used as the prototype of a DSC lowpass filter. 
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The feasibility of the DSC-RSK as lowpass filters in numerical computation is based 
on an important observation: by adjusting the ratio r (r = ■^) of the regularizer, one can 
control the effective wavenumber region and thus, control the resolution of the lowpass filter. 
We plot the frequency response profile of the DSC-RSK lowpass filter in FIG. 1. A distinct 
property of the RSK lowpass filter is that there is no ripples in its passband; thus the 
application of this filter will not impart the low frequency. In contrast, some lowpass filters 
constructed with certain polynomials may have ripples in their passband, which would rule 
out their application to numerical integrations, especially for long time evolution, where the 
solution may be distorted after repeated filtering. 

The maximum effective wavenumber (frequency) range of the DSC-RSK is determined by 
the ratio r, which in turn is determined by computational bandwidth W. The philosophy 
is simple since with the increase of W, the Gaussian window should also be enlarged to 
fully utihze the stencil. The optimal values r at different W can be estimated and they 
characterize the effective wavenumber range. The larger the optimal r is, the wider the 
effective wavenumber range will be. When r value is chosen larger than the optimal value 
for a given W, the accuracy of the resulting filter will be degraded by oscillations. Hence 
we do not recommend the use of r beyond the optimal value unless in some circumstances 
where the filer is utilized to stabilize a long-time integration. 

B. Fourier domain algorithm 

We denote C{xj) — Sa,A{xj) the proposed DSC-RSK lowpass filter in the physical domain 

[0, L] and C(cc;„) its image in the Fourier domain, where A = L/N, Xj = jA and ujn = 
The obvious relation between Ci^j) Ci^n) is 

N-l 

CM = Aj2C{xj)e-'^-''^. (8) 

3=0 

As in the previous work [15], filtering in the Fourier domain is usually accomplished by 
multiplying the given Fourier coefficients u{uJn) = AY^^^q u{xj)e~^'^"'^^ by C{wn) directly, 
i.e., a windowed Fourier transform, 
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Therefore, the modified Fourier coefficients ({u;n)u{uJn) in Eq. (9) are the Fourier coefficients 



of a localized version of the original function u{xj): u''{xj). Hence, errors arising from 
the discontinuity are also localized and the accuracy away from the discontinuity can be 
ensured. In practical applications, one only needs to carry out the filtering and differentiation 
D^u — Ux in one step 



A standard forth-order Runge-Kutta (RK-4) scheme is employed for the time advance scheme 
in this work. We use a TVD switch to adaptively activate the apphcation of the DSC- 
RSK lowpass filter. Whenever the total variation of the approximate solution in the two 
consecutive time steps exceeds a prescribed criteria, the lowpass filter is activated. The 
DSC-RSK lowpass filter is only apphed at the end of each full RK-4 cycle. The reader is 
referred to our earlier work [19,45,46] for the detail of this adaptive filtering. 



It is more difficult to apply the filter in the physical domain than in the Fourier domain. 
On one hand, it is not easy to quantitatively analyze the infiuence of numerical dissipation 
in the physical domain. On the other hand, the convolution operation 



complicates the physical domain application of the filter. However, in some practical appli- 
cations, it turns out to be quite necessary to implement the filter in the physical domain, 
especially for problems whose boundaries require some special treatments. The refiective 
boundary is one of such examples that need to be treated by using grid points outside the 
computational domain. In this work, we propose the method to apply our DSC-RSK lowpass 




(10) 



C. Physical domain algorithm 




(11) 



filter in the physical domain. Due to the fact that the RSK lowpass filter is an interpolatory 
kernel, we consider a two-step procedure. Step one is the prediction of mid-mesh values 
from the nodal values 

-1 1^1 
<i = E «.+.+iC((j + ^)A) + E«^+.C((J - x)A). (12) 

2 j=--w ^ j=l ^ 

The other is the reconstruction of the nodal values from the mid-mesh values 

-1 -1 H/ 1 

= E <,+iC((j + ^)A) + E<,-iC((j - ^)A), (13) 

where, Ui denotes u{xi) and A is the grid spacing. The lowpass filter with a r value smaller 
than the optimal value is used in one of these two steps. The activation of the RSK lowpass 
filter is controlled in the same manner as that of the frequency domain algorithm. 

As described, the DSC-RSK lowpass filter can be applied either in the physical domain 
or in the Fourier domain. In general, filtering in the Fourier domain is much easier and 
more convenient than in the physical domain because no additional summation is required 
to implement the filter in the Fourier spectral method. In contrast, the physical domain 
implementation requires two additional summations. Therefore, we have used the filter in 
the Fourier domain in most of our numerical experiments except for the cases in which the 
refiective boundary condition is involved. We have checked that whenever applicable, results 
obtained from the Fourier domain approach and the physical domain approach are identical. 



D. Post processing filter 

Once the time integration has been completed, a strong filter is usually required for 
post processing the solution so as to make the solution more presentable. In principle, we 
only need to reduce the r value of the DSC-RSK lowpass filter at the last time step of the 
computation. As an alternative, we provide a series of lowpass filters based on the Lagrange 
delta kernel which is given by 

w _ 

5w.(x,x,-)= n j = ,w (14) 
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FIG. 2 shows the magnitude response of the Lagrange delta kernel. In fact, the highpass 
versions of these filters are graphically identical to centered finite difference (FD) schemes. 
In our numerical experiments, we adopt Lagrange-2 or Lagrange-4 as post processing filters. 
Their efficiency is demonstrated by our numerical results. 

III. NUMERICAL EXPERIMENTS 

In this section, we examine the validity and demonstrate the performance of proposed 
windowed (filtered) Fourier spectral scheme for conservation law systems. A large collection 
of standard linear and nonlinear benchmark problems [36], including the linear advection 
equation, ID inviscid Burgers' equation, shock tube (Sod and Lax) problems, ID shock- 
entropy interaction, 2D shock-entropy interaction, 2D shock-vortex interaction, 2D advec- 
tion of an isentropic vortex and flow past a cylinder, are considered in the present study. 
Some of these problems have periodic boundary conditions therefore the windowed Fourier 
pseudospectral method will be applied directly. For problems whose boundary conditions are 
not periodic, we symmetrically double the computational domain in both x and y directions 
to compute the flux derivative. It is noted that the extension is necessary for construct- 
ing periodic problem feasible for the Fourier spectral method. Our filtering procedures, 
as described earlier, can be applied with general boundary conditions, such as Dirichlet or 
Neumann. In the appendix, we list the filter parameter r used for each numerical example. 

A. Scalar conservation law systems 

We begin our numerical experiments with the scalar conservation law system, which is 
given by 

«t + /(«)x = 0, (15) 

where f{u) is a function of u. It is generally believed that spectral methods are not very 
suitable for simple shock profiles and should not be used if one is interested primarily in the 
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shock profiles rather than other detailed structures of the fiow. In this study, we choose scalar 
conservation law systems to illustrate that the proposed global spectral scheme performs well 
for a class of localized shock problems. 

Example 1. We first solve a linear advection equation, which is of the form 



(16) 



Ut -\- Ux — ^ — 1 < X < 1, 

u{x, 0) = Uo{x) periodic, 
where, uq is an initial value 

l{G{x, (3,z-6) + G{x, p,z + 6)+ 4G{x, (3, z)) -0.8 <x< -0.6 

1 -0.4 <x< -0.2 

tio(a;) = "i 1 - |10(x- 0.1)1 0<x<0.2 (17) 

l{F{x, a,a-6) + F{x, a,a + S)+ 4F(a;, a, a)) 0.4 < a; < 0.6 

otherwise. 

Here, functions G and F are defined as 



G{x, P, z) — e 



(18) 

F{x^a^a)= Jm.ax.{l — a^{x — aY^Q)^ 

where a = 0.5, z = -0.7, 5 = 0.005, a 10 and /? = 

The solution of this problem contains a combination of a Gaussian, a square wave, a 
sharp triangle wave and a half ellipse. The exact solution at any time can be easily obtained 
by a translation of the initial solution at speed 1. In FIG. 3 and FIG. 4, we display the 
numerical results for this problem at t=8 with 128 and 256 grid points, respectively. We 
observe that our method works well in both cases. 

Example 2. Secondly, we test the proposed method by considering a moving W-shape 
wave, i.e., a piecewise continuous initial value, for the linear advection equation with an 
initial value [45] 
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uo{x) = < 



4x-l 



-4x + f 



1 





< a; < 0.2 
0.2 < X < 0.4 
0.4 < X < 0.6 
0.6 < a; < 0.8 
otherwise. 



(19) 



This case is very similar to Example 1. It contains the so-called contact discontinuity and 
is quite difficult to solve in hyperbolic conservation laws. In particular, the interaction of 
two lines at x — 0.2 and x — 0.6 is difficult to resolve. Hence, it is a good test for the 
shock-capturing abihty of the present method. We compute the solution up to t—S. The 
results are shown in FIG. 5 and FIG. 6. 

Example 3. Having solved the linear equation with different initial values, we consider 
the solution of the most popular inviscid Burgers' equation, in which f{u) — v? /2 and the 
Riemann type initial value is given by 



u{x, 0) 



(20) 



1 X < 
a; > 0. 

This is a standard benchmark problem in hyperbolic conservation laws and has been consid- 
ered by numerous researchers. The exact solution is a shock wave with a constant velocity 

1 x-St<{) 

u{x, t) = < 

x-St>0, 

where the speed of the shock front is 



(21) 



1 

2' 



(22) 



We should note that this case is a non-periodic boundary problem. We must symmetrically 
double the computational domain in the x-direction when we calculate the derivative of the 
flux. This operation is necessary because it generates the periodic boundary condition so 
that the Fourier spectral method can be used. The numerical results are plotted in FIG. 7 
at time t—2. 
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Example 4- Another Riemann type initial value for the inviscid Burgers' equation with 



a flux oi is given by 



u{x, 0) 



X <0 

1 x>0. 



(23) 



The exact solution of this problem is a rarefaction wave 



u{x, t) 



f < 



< f < 1 
f >1- 



(24) 



Numerical results are plotted in FIG. 8 at t—2 with different numbers of nodal points. 
Particularly, one can see that the end of the rarefaction fan is well resolved. 

Example 5. Finally, we use non- convex flux to test the convergence to the physically 



correct solution. The flux is a non-convex function 



with a Riemann type initial value 



(25) 



u{x, 0) 




(26) 



The reader is referred to [20] for more detailed information about this problem. The re- 
sults are displayed in FIG. 9 at t—O.OA. Again, present spectral solver yields a satisfactory 
resolution. 

B. Euler systems in one dimension 

In this subsection, we perform numerical experiments by using the proposed scheme for 
the ID Euler equation of gas dynamics. In one dimension, the Euler equation takes the form 



Ut + F{U\ = 



(27) 
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with 



U 



( \ 

P 

pu 



pu 

,2 



F{V)^ pu-'^p , (28) 

u{E+p) 

where, p, p and E denote the density, velocity, pressure and total energy per unit mass 
E — p{e + (it^)/2), respectively. Here, e is the specific internal energy. For an ideal gas 
with the constant specific heat ratio (7 = 1.4) considered here, one has e = p/(7 — l)p. We 
consider the following two well-known Riemann problems. 

Example 6. We solve for the solution of two shock tube problems with Sod's and Lax's 
], i.e.. 



initial conditions 

{p,u,p)t=o = < 
for the Sod problem and 



:i,o,i) 



x < 



(0.125,0,0.1), otherwise 



(29) 



t=0 



(30) 



(0.445,0.698,3.528), x<0 
(0.5,0,0.571), otherwise 

for the Lax problem. The numerical results for these two problems are shown in FIG. 10 and 
FIG. 11, respectively. We can see that the present method gives a good resolution in both 
cases. Similar to all the methods depending on the numerical dissipation, the resolution 
for the linear degenerate contact discontinuity is slightly lower than that for the generic 
nonlinear shock wave. If combined with Marten's artificial compression method, present 
spectral solver should have better resolution for the contact discontinuity. 

Example 7. In this example, the interaction of an entropy wave of small amplitude with 
a Mach 3 right-moving shock in a one-dimensional fiow is investigated. The computational 
domain is taken as [0, 9] and the fiow field is initialized with 



(p, u,p)t=o 



(3.85714, 2.629369, 10.33333) x < 0.5 

(^^-eSininx)^ 0, 1.0) x>Q.5 



(31) 



14 



where e and k are the amphtude and the wave number of the entropy wave before the shock. 
This problem is significant due to its relevance to the interaction of shock-turbulence. Spec- 
tral methods, by nature of their high accuracy and negligible dispersive and dissipative 
errors, are ideally suited to the numerical study of turbulence and have already been widely 
applied. Here we would hke to argue the feasibihty of the present spectral method to the 
shock-turbulence interaction. In our test, we vary the wavenumber of the pre-shock wave 
while keeping its amplitude unchanged (e = 0.01). As the wavenumber increases, the prob- 
lem becomes more and more challenging because the amplified high-frequency entropy waves 
after the shock are always mixed with spurious oscillations. As the wavenumber increases, 
it is very difficult to distinguish the high-frequency waves from the spurious oscillations. 
A low-order scheme may dramatically damp the transmitted high frequency waves. Even 
some popular high-order schemes also encounter the difficulty in preserving the amplitudes 
of the entropy waves due to their excessive dissipation with a given mesh size. Therefore, 
a successful shock-capturing method should be able to eliminate Gibbs' oscillation while 
capturing the shock and preserving the high-frequency entropy waves. 

In our test, we let the shock move from x — 0.5 io x — 8.5. For the purpose of comparison 
with the previous results, we only give the results at the interval [4.0,9.0]. Also, in order 
to discharge the transient waves due to the non-numerical initial shock profile, we plot the 
length of the amplified entropy waves in the same manner as that in Ref. [23]. Furthermore, 
we must point out two nontrivial remarks for the following cases with different k. One is 
that we have replotted the final results in a denser grid, otherwise the bounded strip will 
not be fully spanned due to the fact that no enough grid values locate near all the peaks 
/valleys of the wave. Another is that we have used a local filter in order to eliminate the 
oscillations near the shock when we present final results after the computation. 

First, we consider k — 13. In this case, the computational domain is deployed with 513 
grids points. Such a mesh is suitable for the Fourier spectral method. The numerical results 
of the amplitude of the entropy waves are shown in FIG. 12(a). It can be seen that the shock 
and the generated entropy waves are ideally captured. The entropy waves fully span the 
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strip bounded by two solid lines at ± 0.08690716, which show the amplitude of the amplified 
entropy waves predicted by the linear analysis. Since post-shock waves are monochromic, a 
simple calculation indicates the point per wavelength (PPW) value of these waves is 7.125. 
Such value is larger than the spectral resolution (PPW=2.0) but it is much smaller than 
that obtained by most shock-capturing schemes. The transmitted waves suffer the similar 
smearing as the contact waves. The shock wave, on the contrary, is self-repairable due to 
its compression nature. This observation verified that the numerical dissipation needed for 
shock-capturing will degrade the resolution of the spectral method. It also verified that this 
degradation can be isolated and alleviated by the DSC-RSK filter. 

We next double the wave number k up to 26. A mesh of 513 grid points is not enough for 
simultaneous shock capturing and high-frequency wave resolving. Therefore, we increase the 
mesh size to iV = 1025 for this computation. The results for this case are presented in FIG. 
12(b). It is seen that no obvious deterioration happens to the results and the resolution in 
the post-shock waves is still satisfactory. The fullness of the wave is comparable with that 
of «; = 13. But the first two post-shock waves have been slightly polluted. 

When the wave number k, is increased to 39 and 52, we use 2049 grid points to resolve 
such high-frequency waves. FIG. 12(c) and FIG. 12(d) show the results of k = 39 and 
K — 52, respectively. The same PPW value 7.125 is maintained in this high-frequency case. 
The good resolution of the present method in simulating this ID shock entropy interaction 
indicates that our scheme is capable of and efficient in distinguishing the high-frequency 
entropy waves from spurious oscillations. 

Example 8. Results are now shown for the problem of Shu and Osher which also de- 
scribes the interaction of an entropy sine wave with a Mach 3 right-moving shock. The 
computational domain is taken as [—1, 1] and the fiow field is initialized with 

f (3.85714, 2.629369, 10.33333) x < -0.8 

{p,u,p)t=o=\ (32) 
I (1.0 + esin(«;7rx), 0, 1.0) x > -0.8, 

where e=0.2 and K—b. FIG. 13(a) and FIG. 13(b) show the solution of the density with 128 
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and 256 cells, respectively. The 'Exact' solution here is the solution computed by the third- 
order ENO with 1200 cells. As evident in the solution, the complicated flow fleld behind 
the shock is well resolved and the shock remains sharp. However, small oscillations near the 
shock still exist due to the fact that we do not attempt to locate this shock intentionally. If 
we predict the position of the shock and apply a local fllter near the shock at the final time 
step, it is believed that more satisfactory results can be obtained. 

C. Euler systems in two dimensions 

In two dimensions, the Euler equation for gas dynamics in a vector notation takes the 
conservation form 



Ut + F{U), + G{U)y = Q 



(33) 



with 



U 



( \ 
P 

pu 
pv 



I 



■ F{U) 



pu 
pu^ + p 

puv 
u{E + p) 



; G{U) 



pv 
puv 
pv^ +p 
v{E + p) 



(34) 



p^{j-l){E--p{u' + v')), 



(35) 



where {u, v) is the 2D fluid velocity and 7 = 1.4 is used in all computations. 

Example 9. As the flrst example, we use our spectral scheme to study the interaction 
between a normal shock and a weak entropy wave which makes an angle 9 G (0, 7r/2) against 
the X-axis. If ^ = 0, we have essentially the ID problem (see example 7). The initial 
conditions are deflned as follows: the right state of the shock is given as {pr,Ur,Vr,Pr)—ii, 
0, 0, 1) with a given shock Mach number Mg = 3. We add a small entropy wave to the flow 
on the right of the shock which is equivalent to changing only the density of the flow on the 
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right of the shock. In our problem we change the density pr in the right state of the shock 
by multiplying 

"= (36) 

„ „—{e/pr)sm.{K{xcosd+ysin&)) 

Pe — 

with e and k being the amplitude and wavenumber, respectively. In order to carry out a 
long time integration and to enforce periodic boundary condition in y-direction, the com- 
putational domain is taken to be [0,9] x [0, -;^^\- We initially position the normal shock 
at X = 0.5 and allow it to move up to x = 8.5. In our simulation, we adopt e = 0.1, 
K — lb, and 6 = 30°. The performance is measured by drawing the maximum amplitudes 
of the amplified entropy waves in the y-direction for all fixed grid values x e [7.4, 8.4], and 
comparing them with the amplitude predicted by the linear analysis, which is 0.08744786. 

We use 513 grid points in the x-direction and 32 grid points in the y-direction in our test. 
In order to avoid peak/ valley losses of the entropy waves due to insufficient grid points near 
these positions, we interpolate our final results to a denser grid when plotting them. FIG. 14 
displays the amplitudes of the amplified entropy waves. Obviously, the compressed entropy 
waves are well resolved in the post-shock regime even though small oscillations remain near 
the location of the shock. It is believed that such a trivial flaw is acceptable. 

This example has further demonstrated the capability of our spectral scheme for captur- 
ing small scale waves in the presence of a shock. 

E'x ample 10. This model problem describes the interaction between a stationary shock 
and a vortex. It has many potential apphcations and has been treated by many researchers 
using various techniques. Early numerical studies in this area using spectral methods focused 
primarily on shock-fltting techniques. However, for the cases in which a vortex causes a shock 
to deform at the point of bifurcation, it is generaUy difficult for a shock-fitting method to 
apply. This problem was successfully simulated in Refs. [6,8] by using the Chebyshev spectral 
method. In this paper, we also adopt this model to access the capability of the proposed 
spectral shock-capturing method. This problem is defined as follows. The initialization is 
imposed on the domain [0, 2] x [0, 1] with a stationary normal shock x — 0.5 and a Mach 
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1.1 flow at the inlet. The right state of the shock is given as {pr, Ur, Vr,Pr) — (1, 1-1 0> 
A vortex is generated by a perturbation to the velocity («, v), temperature T and entropy 
S which is centered at (xc, Uc) = (0.25, 0.5) and in the forms 

ii' = ere°(i-^')sin^ (37) 
^;' = -ere"(^-"')cos^ (38) 

r = - ^^ ' (39) 

4q;7 

5' = (40) 

where, r = ^, r = \j{x — + {y — Z/c)^ and 9 = tan~^{{y — yc)/{x — Xc))- Here e denotes 
the strength of the vortex, a is the decay rate of the vortex. Note that we can obtain 
perturbations in p and p according to the relations of T — p/p and 5" = In^. In our 
numerical experiment, we choose e = 0.3, Tc = 0.05 and a = 0.204. 

A reflective boundary is imposed on both the upper and lower boundaries. We solve 
the Euler equations on a 257 x 129 Cartesian grid, which is uniform in y-direction and 
refined in x-direction around the shock using the Roberts transformation [2]. Once again, 
we emphasize that the DSC-RSK lowpass filter is used in the physical domain directly for 
this problem due to the presence of the refiective boundaries. FIGs. 15(a-e) show a time 
evolution of shock- vortex interaction in terms of the pressure contours. We can see that even 
at t—0.8 the solution (FIG. 15(e)) is almost free of oscillations and the fine scale features 
are captured very well. Particularly, the shock bifurcation reflected on the top boundary is 
shown clearly. 

In order to explore the potential of our algorithm further, we increase the flow Mach 
number up to 3.0 and the strength of the vortex e up to 0.6. The results at t—0.2 are 
shown in FIG. 16. Deformation and bifurcation of the shock are observed clearly. These 
results indicate that the applicability of our spectral shock-capturing scheme is not limited 
to low Mach number and low vortex strength. Flows of high Mach number and high vortex 
strength can also be successfully simulated. 

Example 11. In this case, we investigate the advection of an isentropic vortex in a free 
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stream. As the exact solution of this problem is available, it is a good test for examining 
the accuracy and stability of shock capturing schemes. 

Consider a mean flow of (poo, '"oo, i^oo,Poo, ^oo) = (1, 1, 1, 1, 1) with a periodic boundary 
condition in both directions. At to, the flow is perturbed by an isentropic vortex {u',v',T') 
centered at {xo,yo), having the forms of 

u'^-^{y-yo)e'^'-''\ (41) 

v'^ ^{x - xo)e^^'-''\ (42) 
j..^_(7^_l)A!g2,(i-.^)_ (43) 



Here, s — \j{x — XqY + (y — VoY is the distance to the vortex center; A is the strength of 
the vortex and is a parameter determining the gradient of the solution. In our test, we 
choose them as A = 5 and r] — 1 unless we further specify. Note that for an isentropic flow, 
relations p — p^ and T — p/ p are valid. Therefore, the perturbed p is required to be 



(7 - 1)A ^2r,(l-.^) 



(44) 



The computational domain is taken to be [0, 10] x [0, 10] with the center of the vor- 
tex being initially located at (xo,yo) = (5,5), the geometrical center of the computational 
domain. 

Since there is no presence of shock in this problem, spectral method on its own can 
already provide excellent results if the time integration period is small enough. As such, the 
lowpass filter is not needed in an initial time period. However, as time progresses, errors 
would accumulate rapidly and the computation could become unstable if the lowpass filter 
is not activated to efficiently control the dramatically nonlinear growth of the errors. In our 
numerical experiment which is not shown here, the computation blows up at around t=13 
if no filter is used. 

Two experiments are conducted in this study. One is to examine the accuracy of our 
spectral code and the other is to investigate the stability of our scheme for a long time 
integration. For the first experiment, we compute the solution of p at the time t—2 using 
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three sets of meshes {N — 32, 64, 128). To make the spatial error dominant, we optimize the 
CFL number to be 0.01. 

Two error measures, Li and L2, have been used. To be consistent with the hterature 
[11], two errors used in this problem are defined as 

N N 

(^^ + ~.=o,=o 
1 



~ (N + 1)2 I fiJ - fi,j I (45) 



N N 



EEI/m-/mP, (46) 



where / is the numerical result and / the exact solution. Note that they are not the standard 
definitions. The errors for the density with respect to the exact solution are listed in Table 
I in which we also hst the results of some other schemes for a comparison. As the table 
highlights, the accuracy of the Fourier spectral method is comparative to that of the DSC 
method [48], while both Fourier spectral and DSC methods are much more accurate than 
the fourth-order accurate, conservative central scheme (C4). 

Our next numerical experiment concerns the performance of our spectral code for a long 
time integration. Here, we compute the solution of p at the time t—100, 200, 400, 600, 800 
and 1000, with the grid — 64. Note that we choose the gradient of the solution rj — 0.5 
in this case. We list the Li and L2 errors for the density in Table II. It can be seen that our 
spectral code is highly accurate even after such a long time evolution. The results indicate 
that our spectral scheme together with the DSC lowpass filter is reliable for the numerical 
approximation of conservation laws. 

Example 12. Finally, we consider the problem of a supersonic fiow past a cylinder [23]. 
This example is used to mainly examine the ability of our DSC filter spectral scheme in 
handling a non-rectangular domain. 

The physical domain is on the x — y plane, and the computational domain is chosen to 
be [0, 1] X [0, 1] on the ^ — 77 plane. The mapping between the computational domain and 
the physical domain is 
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X = -{R, - {R, - mcos{9{2r] - 1)) 

(47) 

y = {Ry- {Ry - msm{e{27j - 1)) 

where we take R^ — Ry — 6 and ^ = ff ■ 

The flow fleld is initiahzed with a Mach 3 shock moving from the inlet, i.e., = 0, toward 
the cyhnder. The state in front of the shock is given as {pr, Ur, Vr,Pr) — (1-4, 1.0, 0, 1.0). We 
impose a reflective boundary condition on the surface of the cyhnder, i.e., ^ = 1, and the 
outflow boundary condition at ?7 = and f] = 1. In our simulation, a uniform mesh of 65 x 129 
in the computation domain is used. FIG. 17(a) depicts a diagram of the mesh (drawing every 
other grid line) in the physical domain. Like the case on a rectangular domain, we double 
the computational domain in both ^ and 7] directions when we approximate the flrst order 
derivative of the flux using the Fourier spectral method. The numerical result given in terms 
of pressure contours (obtained at t—A.5) is shown in FIG. 17(b). In summary, our spectral 
scheme can be easily applied to the nonrect angular domain, as long as the domain can be 
smoothly transformed to the rectangular one. 



IV. CONCLUSION REMARKS 

In this work, we introduce the use of discrete singular convolution (DSC) fllters in the 
Fourier pseudospectral method (FPM) for the simulation of hyperbolic conservation law 
systems. The fast Fourier transform (FFT) is used as the basic scheme, while DSC low- 
pass filters arc used to suppress unphysical oscillations. The DSC filters are implemented 
in either the Fourier domain or the physical domain, depending on the physical problem 
under consideration. The Fourier domain implementation is easy and cost efficient, whereas 
the physical domain implementation allows more fiexibility to handle some special bound- 
ary conditions. A sensor technique developed in the conjugate filter oscillation reduction 
(CFOR) scheme [19,45,46] is adopted in the present work to appropriately activate the DSC 
filter. The fourth-order Runge-Kutta scheme is utilized for the time integration. 

Extensive numerical experiments are considered to validate the proposed approach and 
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to demonstrate its usefulness. Excellent numerical results are obtained for all problems ex- 
amined. The DSC filters turn out to be very effective in removing the Gibbs oscillations 
produced by the spectral approximation while maintaining a high accuracy. Moreover, it 
is also highlighted that these filters can dramatically improve the stability of the spectral 
method for the long time integration. The main features of proposed method are the follows: 

(1) The basic scheme, the FPM is of spectral accuracy, which is important for problems that 
require not only capturing shock waves, but also resolving fine flow structures. 

(2) The use of FFT endows the method high efficience, which is desirable for large scale 
problems in science and engineering. 

(3) Numerical dissipation of the DSC filters is adjustable for a given problem, which endows 
the proposed method the ability to handle problems that are sensitive to numerical dissipa- 
tion. The interaction of turbulence and shock is arguably a problem of this kind. 

(4) The Fourier domain algorithm can be regarded as a windowed Fourier pseudospectral 
method for shock-capturing. The alternative use of the DSC filter in the physical domain 
makes the proposed method feasible to problems requires special boundary conditions, e.g., 
reflective boundary condition. 

(5) The proposed scheme by-passes finding the characteristics which is required in many 
other shock-capturing schemes. 

(6) The use of the method in higher dimensions is straitforward. 

It is note that the state of art shock-capturing schemes, such as WENO and central 
schemes, do not depend on any adjustable parameter. To make the proposed scheme robust 
and user friendly, future work will focus on developing sophisticated sensors so that the filter 
parameter will be chosen automatically for each given problem. We believe that this can 
be accomplished by analyzing the Fourier power spectrum distribution of the approximate 
solution at selected time steps. 
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FIGURES 
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Wavenumber 

FIG. 1. Magnitude response of the DSC-RSK filter, r = 1.0,1.5,2.0,2.5,3.2,5.0 from the left 



to right. 
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FIG. 2. Magnitude response of the Lagrange filter. Lagrange-2, Lagrange-4, Lagrange-8 and 
Lagrange- 16 from the left to right. 
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FIG. 5. Linear advection equation. t=8, At=0.001, 128 grid points. 
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FIG. 6. Linear advection equation. t=8, At=0.001, 256 grid points. 
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(a) (b) 
FIG. 7. Inviscid Burgers' equation {t=2, A t=0.005). (a) 129 grid points; (b) 257 grid points. 
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(a) (b) 
FIG. 9. Inviscid Burgers' equation with non-convex flux (t=0.04, At=0.0005). (a) 65 grid 
points; (b) 129 grid points. 





(a) (b) 
FIG. 10. Sod's problem (t=2, At=0.02, 129 grid points), (a) Density; (b) Pressure. 
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(a) (b) 
FIG. 11. Lax's problem {t=1.5, At=0.02, 129 grid points), (a) Density; (b) Pressure. 



34 




-0.1 



(a) K = 13, AT = 513 




4 5 6 7 8 9 

(c) K = 39,N = 2049 



35 



0.1 
0.05 


-0.05 
-0.1 




(d) K = 52, TV = 2049 
FIG. 12. ID shock entropy wave interaction. Amplitude of entropy waves. 
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(a) (b) 
FIG. 13. Shu-Osher problem (i=0.47). (a) 128 cells; (b) 256 cells. 
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FIG. 14. 2D shock entropy wave interaction {N = 513). Amplitude of entropy waves. 
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(a) t=Om, 30 contours (b) i=0.20, 30 contours 




(c) f=0.35, 30 contours (d) t=0.60, 90 contours from 
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(e) i=0.80, 30 contours from 1.01 to 1.29 
FIG. 15. 2D shock vortex interaction (Mach=l.l, e = 0.3, CFL=0.5). Pressure profiles. 
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TABLES 



mesh 


error 


FFT 


C4 


DSC 


N=32 


Li 


5.58E-5 


1.96E-3 


1.27E-4 




L2 


1.27E-4 


5.67E-3 


4.07E-4 


iV=64 


Li 


2.33E-8 


1.32E-4 


8.18E-8 




L2 


7.94E-8 


4.39E-4 


4.07E-7 


A^=128 


Li 


4.01E-11 


8.90E-6 


3.99E-11 




L2 


5.09E-10 


2.96E-5 


5.09E-10 



TABLE I. Advective 2D isentropic vortex. Li and L2 errors for the density at t=2. The CFL 
number is 0.01 for all schemes. C4 denotes the fourth-order accurate, conservative centered scheme. 



mesh 


error 


t=100 


t=200 


t=400 


t=600 


t=800 


t=1000 




Li 


4.63E-7 


5.05E-7 


l.OOE-6 


1.44E-6 


2.07E-6 


2.77E-6 


A=64 


















L2 


1.99E-6 


1.23E-6 


2.90E-6 


4.42E-6 


6.02E-6 


7.59E-6 



TABLE IL Advective 2D isentropic vortex. Li and L2 errors for the density at different times 
with iV=64 (CFL=0.5, ?7=0.5). 
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APPENDIX 

The optimal values (r) for DSC-RSK lowpass filters in our numerical experiments. 



No. of Example 


case 


r 


1 


iV=128 


0.6 




N=256 


0.8 


2 


N=128 


0.6 




iV=256 


0.8 


3 


A^=129 


0.7 




N=257 


0.8 


4 


N=129 


0.6 




N=257 


0.6 


5 


N=65 


0.8 




N=129 


0.8 


6 


Lax 


0.95 




Sod 


1.1 




K=13 


2.0 


7 


K=26 


2.0 




K=S9 


2.1 




K=52 


2.1 


8 


N=129 


2.0 




A/"— 257 


2.1 


9 




2.1 


10 




2.8 


11 


r/=1.0 


3.2 




?7=0.5 


2.8 


12 




1.2 
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